Dipolar fermions in a two-dimensional lattice at non-zero temperature 
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We examine density ordered and superfluid phases of fermionic dipoles in a two-dimensional 
square lattice at non-zero temperature. The critical temperature of the density ordered phases 
is determined and is shown to be proportional to the coupling strength for strong coupling. We 
calculate the superfluid fraction and demonstrate that the Berezinskii-Kosterlitz-Thouless transition 
temperature of the superfluid phase is proportional to the hopping matrix element in the strong 
coupling limit. We finally analyze the effects of an external harmonic trapping potential. 



INTRODUCTION 



An increasing number of experimental groups are trap- 
ping and cooling atoms or molecules with a permanent 
magnetic or electric dipole moment. Bose-Einstein con- 
densates of 52 Cr atoms [TJ [2] and of 164 Dy atoms [3 
with large magnetic dipole moments have been realized. 
Fermionic gases of 40 K 87 Rb [1] and 23 Na 6 Li [5] molecules 
with an electric dipole moment have been created, and 
the first steps toward the formation of fermionic 23 Na 40 K 
molecules have been reported [BJ. Also, experimental 
progress toward realizing dipolar molecules in an optical 
lattice have recently been presented [7J. The anisotropy 
of the dipole interaction results in many intriguing ef- 
fects. In a two-dimensional (2D) lattice, the existence of 
density ordered phases with a complicated unit cell [8], 
liquid crystal phases [5] , and a supersolid phase [TO] have 
been predicted when the dipole moments are perpendic- 
ular to the lattice plane. Tilting the dipoles toward the 
lattice plane leads to density order with different symme- 
try, superfluidity and bond-solid order at zero tempera- 
ture [TTHT3"] . When a trapping potential is present, these 
phases were shown to coexist, forming ring and island 
structures [T2"] . 

In this paper, we examine fermionic dipoles in a 2D 
square lattice including the presence of a harmonic trap- 
ping potential. Focus is on the effects of a non-zero tem- 
perature and the melting of density ordered and super- 
fluid phases. We determine the critical temperature for 
the density ordered phases and find that it is propor- 
tional to the interaction strength in the strong coupling 
regime. For the superfluid phase, we calculate the su- 
perfluid fraction and the Berezinskii-Kosterlitz-Thouless 
(BKT) transition temperature, which is proportional to 
the hopping matrix element in the strong coupling limit. 
We analyze the effects of an external trapping potential 
showing that for experimentally realistic systems, the or- 
dered phases exist in the center of the trap with melting 
temperatures close to that which can be obtained from a 
local density approximation. 



MODEL 

We consider fermionic dipoles of mass m and dipole 
moment d moving in a 2D square lattice with lattice 
constant a. The dipole moment is aligned by an external 
field to form an angle dp with respect to the z axis which 
is perpendicular to the lattice plane and an angle <pp 
with respect to a lattice vector chosen as the x axis. The 
Hamiltonian is H = -ffkin + V where 

-ffkin = ~*y"l (c\cj + h.c^j + ^ ( ^ m ^ 2 r 2 - fij hi 

(1) 

and 

V = lY,V D (T ij )n i fij (2) 

where Yi denotes the position of lattice site i and r.y = 
Yi — Yj, di is the annihilation operator that removes a 
dipole at site i, and hi = c|cj is the number operator. 
The chemical potential is /i and t is the hopping matrix 
element between nearest neighbors (ij). We include the 
effects of a harmonic potential with trapping frequency 
u) exactly in our analysis. The interaction between two 
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FIG. 1: (color on-line) Dipoles move in a 2D square lattice 
with lattice constant a. They are aligned forming an angle 
Op with the z axis perpendicular to the lattice plane, and 
the azimuthal angle 4>p with the x axis which is parallel to a 
lattice vector. 
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dipoles separated by r is given by 



Vd{v) 
D 2 



D 2 



1-3 cos 2 



[1 - 3 cos 2 Op 



i sin 2 0p] 



(3) 



with D 2 = d 2 /47re for electric dipoles and # r d the angle 
between d and r = r(cos0,sin0,O), see Fig.[l] We define 
g = D 2 /a 3 as a measure of the interaction strength. 

The anisotropy of the dipolar interaction ^ with at- 
tractive and repulsive regions gives rise to both density 
ordered and superfluid phases [Til H2] • We capture the 
existence of these competing phases using mean-field the- 
ory including the Hartree terms and the pairing terms, 
which we expect to be reasonably accurate due to the 
long range nature of the interaction. The mean-field 
Hamiltonian is diagonalized by solving the Bogoliubov-de 
Gennes equations [12] 



E 



(4) 



where A,.,- = Vd(iYj )(cj-Cj) and 



Lij = -tS {ij) + Vb(rifc)(«fc) + yw 2 r- - fi)S ir (5) 



Here <5^ and 6^ are the Kronecker delta functions con- 
necting on-site and nearest neighbor sites, respectively. 
Self-consistency is obtained iteratively through the usual 
relations (h t ) = T,e„>o [(1 - fr,)\vj,\ 2 + f v \u),\ 2 ] and 



(ciCj) = Eb^oK^^ 1 _ fv) + v n u vfv\i with fv = 
[exp(£L/T) + 1] _1 the Fermi function for the temper- 
ature T. We use units where fcg = H = 1. To analyze the 
melting of the superfluid phase, we shall use the frame- 
work of BKT theory. 



STRIPE MELTING AT HALF FILLING 

Wc first analyze the case of no trapping potential and 
half filling, N/N L = 1/2, with N L the number of lat- 
tice sites and N — Ei(^«) the total number of particles. 
When the dipoles are perpendicular to the lattice, it fol- 
lows from the perfect nesting of the Fermi surface that a 
phase with checkerboard density order persists down to 
g/t for T = 8\. In the limit of strong interac- 
tion j/t > 1 where the kinetic energy can be neglected 
and the problem becomes classical, it was shown that 
the checkerboard phase is replaced by a striped phase 
when the dipoles are tilted at a sufficiently large angle 
P [12]. We now examine the melting of these density 
ordered phases at a non-zero temperature. The melting 
is in the Ising universality class due to the discreteness 
of the lattice, and we therefore expect mean-field theory 



to yield a qualitatively correct value for the transition 
temperature. 

For the case of stripes along the x direction, we 
express the density as (h % ) = 1/2 [l + M(-l) Vi ' a ] 
with M the order parameter. The correspond- 
ing mean-field Hamiltonian can be written as H = 



E 



fc„>0 



-Bik7i k 7ik + ^2k72k72k 



with the single particle 



energies 



Ei k = £k + V (2i cos k y a) 2 + [V D (0, ir/a)M /2] 2 (6) 

where £k = — 2tcosk x a— /x+Vd(0,0)/2. We have defined 
the Fourier transform Vb(k) = Ej ex P( — *k ' r i)^D( r i)- 
The energy i?2k is given by ([6| with a minus-sign in front 
of the square root. The self-consistency equation reads 



Mr 



Vb(0,7r/q)(/ lk -/ 2k ) 
Nl k y >o J(2tcosk y ) 2 + [V D (0,ir/a)M/2] 2 



(7) 



where the sum is over half the first Brillouin zone with 
ky> 0. In the limit of strong interaction g/t 3> 1, Eq. 
(m) yields 



~Vb(0,7r/o). 



(8) 



When the dipoles are aligned in the lattice plane with 
{9 pi (f> P ) = (tt/2,0), Eq. H gives Tf ps 1.27<?. A sim- 
ilar analysis for the checkerboard phase yields T£ h — 
—Vd(it /a,ir /a)/4 in the strong coupling limit, which 
gives T c cb 0.66 5 for 6 P = [8]. 

Figure [2] shows the critical temperature as a func- 
tion of the interaction strength for the checkerboard 
phase with Op — and for the striped phase with 
(0p,4>p) = (tt/2,0). The o's and x's are numerical re- 
sults for the stripe and checkerboard phases respectively, 
obtained from solving Q , and the lines the analytical re- 
sults for the strong coupling limit discussed above. Finite 
size effects of the system arc eliminated by neglecting the 
high temperature tail of the order parameter. For exam- 
ple, for the lower right inset in Fig. [2] the elimination 
of the high temperature tail gives the critical temper- 
ature T c cb /i = 0.4. We see that the numerical results 
agree well with the strong coupling results for g/t ^> 1 
whereas the critical temperature becomes exponentially 
suppressed in the weak coupling limit. Note that the 
critical temperature of the striped phase is almost twice 
that of the checkerboard phase, which makes it easier 
to observe experimentally. The upper left inset shows 
how the striped order parameter M decreases with T for 
(9p,(j)p) = (tt/2,0) and g/t — 3.3, and the lower right 
inset shows the checkerboard order parameter M as a 
function of T for Op = and g/t = 1. 

Figure [3] shows the critical temperature of the striped 
and the checkerboard phase as a function of (6p,4>p) 
in the strong coupling regime. It is obtained from 
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FIG. 2: (color on-line) The critical temperature of the striped 
phase for (8p,<fip) = (tv/2,0) (o's) and of the checkerboard 
phase for 6p = ( x 's) as a function of coupling strength for 
half filling obtained from a numerical calculation on 30 x 30 
lattice sites. The dashed lines give the strong coupling results 
Tf = -Vd(0,7t/o)/4 with (8p,<j>p) = (tt/2,0) and T c cb = 
— Vb(7r/a, 7r/a)/4 with Bp = 0. The upper left inset shows 
the striped order parameter M for g/t = 3.3asa function of 
T for (9p,4>p) = (7r/2,0). The lower right inset shows the 
checkerboard order parameter M for g/t — 1 as a function of 
T for 9 P = 0. 

max[— y(0,7r/a)/4, — V"(7r/a,7r/a)/4]. For most orien- 
tations of the dipoles, the critical temperature of the 
striped phase exceeds that of the checkerboard phase. 
We note that the upper left corner in the phase-diagram 
shows a negative critical temperature which indicates 
that none of the two phases we explore are stable in this 
region. 

STRIPE AND SUPERFLUID MELTING AT ONE 
THIRD FILLING 

For smaller filling fractions, the system can be in a 
superfluid state with p-wave symmetry for large enough 
9p [TTJ[T2]. This leads to a competition between density 
and superfluid order in analogy with dipoles moving in 
a 2D plane without a lattice [TH H5] . As an example, 
we now consider the melting of the superfluid and the 
striped phase for the filling fraction N/Nl = 1/3 and 
(6p,(pp) — (7r/2,0). For these parameters, mean-field 
theory predicts the system to be superfluid for g/t < 1.15 
and to exhibit stripe order for g/t > 1.15 at T — [12]. 

For the 2D system considered here, the melting of the 



FIG. 3: (color on-line) The critical temperature in units og 
g of the striped and checkerboard phases as a function of 
the dipole orientation (0p,(j>p) for half filling. A dashed line 
marks the boundary between the stripe and checkerboard 
phases, and the region with no density order is bounded by 
another dashed line. 

superfluid phase is of the BKT type with a transition 
temperature determined by the phase stiffness of the or- 
der parameter [151 E3 • The phase stiffness J x associated 
with a phase twist of the superfluid order parameter in 
the x direction is determined from the energy cost 

F e -^y£<5e 2 . (9) 

i 

Here, Fq is the free energy when the phase of the order 
parameter varies by SO between neighboring sites in the 
x direction and Fq is the free energy when there is no 
phase twist [T8]. Associated with the phase twist, we 
define the superfluid fraction p SjX by writing 

F e - F = ^ Ps , x m*v s 2 = ^-tp s , x 5e 2 , (10) 

where v s — SO /2m* a is the superfluid velocity of the 
Cooper pairs with mass 2m*. The effective mass for the 
dispersion — 2£(cos k x a + cos k y a) is m* = l/2ta 2 . Note 
that the superfluid fraction is dimcnsionlcss. Similar ex- 
pressions hold for the phase stiffness J y and the super- 
fluid fraction p s ^ y for the y direction. 

A linear phase twist along the x direction is equiva- 
lent to acting on the Hamiltonian with the unitary gauge 
transformation 

jy e _ e -iSej2iXi/ajj e iSej2i&i/a 
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where xi is the x-coordinate of particle I [15] . We 
have 5® = 289 since the superfluid order parameter in- 
volves two particles so that the gauge transformation 
gives Aij — » Aij exp[i(xi + Xj)S9/a\. The gauge trans- 
formation only affects -ffkm by introducing a phase factor- 



ies ' c, 



£ e ±i<50 cjci± ea . on the hopping terms connect- 



ing neighboring sites in the x direction. Here, e x denotes 
one lattice step in the x direction. Since we only need the 
energy cost to lowest order in the phase twist to deter- 
mine J from Eq. it is sufficient to use perturbation 
theory in 59. Expanding to second order in 69, we obtain 
H e = H + J + f with 



J = -i69t 

t='- 



c]c 



(12) 



1 + (3m* n)^ 1 (2tt)~ 3 J d 3 kdEfkk 2 for a single component 
superfluid |22) . 

From the phase stiffness, we can extract the transition 
temperature as Tbkt = kJ/2 jTHl [T7] where we have 
taken the average J = (J x + J y )/2 to account for the 
anisotropy of the p-wave pairing. Equations ([£])-( f ) give 
J = Np s t/2NL with p a = (p s ^x + Ps,y)/2, and we finally 
obtain 



7T N 7T n s 



(17) 



with the superfluid density defined as n s = Np s /NLa 2 . 

In Fig. |4| we plot Tbkt as a function of the coupling 
strength obtained from Eq. (17). For comparison, we plot 



the mean- field superfluid transition temperature T* . We 
also plot the critical temperature for the stripe phase 
which is the ground state for g/t > 1.15. For weak cou- 



Since the unitary transformation conserves particle num- 
ber, we can take Fq — F = fie — where f2 = F — p,N 
with N the total number of particles [5D]. The linked 
cluster expansion gives [5T] 



r> e - Cl = (f) - | (J 2 ) 



(13) 



where (. . .) denotes the thermal average with respect to 
the untwisted Hamiltonian and we have used that there 
is no current in the untwisted case, i.e. (J) = 0. Mean- 
field theory gives after some lengthy but straightforward 
algebra 



(f) 



t 



•i# e » +<«'-«-)/„ (14) 



7} . I 



and 



"rj "a u r) u a 



ij r}Ot k,l — — l 

X /„(1 - fa) - <«jjtt^*»-«i+' e -V,(l - /«)] • (15) 



The sums in Eqs. (|14[)-( 15 ) are taken over positive as 



well as negative energies, and we have made use of the 
duality (u^, v v ,E v ) f-> (v*, u* , — E v ) of the Bogoliubov-de 
Gcnnes equations. 

When there is no trap, the Bogoliubov-de Gennes 



equations are straightforward to solve and Eqs. ( 10 ) , ( 14 ) , 



and (15) yield 
N ^ 



Ps.x 



2t 

n k cosk x a - — / k (l - / k ) sin 2 k x a 



(16) 

Here are the BCS quasiparticle energies for the p- 
wave paired state, and n k = w k /k + w k (l — /k)- In 
the continuum limit a — > keeping the density N/N^a 2 
constant, this reduces to the usual expression p SfX = 
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FIG. 4: (color on-line) The critical temperature for the su- 
perfluid phase (*'s) and the striped phase (o's) for (9p,4>p) = 
(tt/2, 0) as a function of coupling strength for one third filling 
obtained from a numerical calculation on a 27 x 27 lattice site. 
The x 's give the mean-field superfluid transition temperature 
T* . For illustrative purposes, we plot the critical temperature 
of the superfluid phase even for g/t > 1.15, where stripe order 
suppresses superfluidity. 

pling, the Tbkt approaches T* as expected [23], whereas 
it is significantly lower for stronger coupling. For strong 
coupling, it follows from Eq. ( 17) that the critical temper- 
ature will saturate at Tbkt 
results yield Tbkt 



t. Indeed, the numerical 
0.12i for g/t > 1 as can be seen 
from Fig. |4j Note however that stripe order sets in for 
g/t > 1.15 which suppresses the superfluid order. Like 
the case for half filling, we have T^* ~ g for the critical 
temperature for the striped phase, which is a higher tem- 
perature than the superfluid transition temperature. It 
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is interesting that both critical temperatures, Tbkt ~ t 
and Tj?* ~ g, can be much higher than that of the antifer- 
romagnetic phase for atoms in a 3D lattice, which scales 
as T N ~ t 2 /U in the strong coupling limit with U ^> t 
the on-site interaction [23J [35] . 

In Fig. [5] we plot the superfluid fraction and the near- 
est neighbor order parameter as a function of T for var- 
ious coupling strengths. As usual for a 2D system, the 
superfluid fraction is discontinuous at the critical tem- 
perature. Contrary to a translationally invariant system, 
the superfluid fraction is less than 1, even for T = [2*5] . 
In the inset, we plot the superfluid fraction and the near- 
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FIG. 5: (color on-line) The superfluid fraction p 3 and the 
nearest neighbor pairing \ {ci+ en .Ci)\ as a function of T for var- 
ious coupling strengths. \{ci+ ew Ci}\: Pink x's for g/t = 0.7, 
pink *'s for g/t = 0.8, and pink A's for g/t — 1.5. p s - Blue V's 
for g/t = 0.7, blue o's for p/t = 0.8, and blue o's for g/t = 1.5. 
The numerical calculations are performed on a 27 x 27 lattice 
with one third filling. Inset: The nearest neighbor pairing 
|(ci+ ex Ci}| (green o's) and the superfluid fractions p„ tX (pink 
x's) and p s . y (red *'s) as a function of g for T = 0. 



est neighbor pairing as a function of coupling strength 
for T = 0. We see that ps x 7^ Ps,j/, which follows from 
the anisotropy of the p-wave paring. Note that the super- 
fluid fraction behaves very differently from the pairing as 
a function of the coupling strength 27J. 

We expect correlation effects to decrease the transi- 
tion temperatures of the ordered phases from what is 
predicted in the present paper. Even so, we believe that 
our results are qualitatively correct due to the long range 
nature of the interaction. This includes the scaling of 
Tbkt j Tp, and T^ b for strong coupling. Our results 
therefore present a useful first analysis of the order phases 
of fermionic dipoles in a lattice at non-zero temperature. 



TRAPPED SYSTEM 

The harmonic trapping potential is always present in 
atomic gas experiments. For T = 0, this leads to the co- 
existence of superfluid and density ordered phases form- 
ing ring and island structures |12) . We now investigate 
these effects at a non-zero temperature. 

Figure [6] (top) shows the density and the checkerboard 
order parameter as a function of temperature for the 
dipoles aligned perpendicularly to the lattice plane with 
(9p, 4>p) = (0, 0). We have chosen ui — cua^m/t = 0.24, 
g/t = 1, and fi/t = 4.23 for the numerical calculations, 
giving N — 207 — 210 dipoles trapped and an average fill- 
ing fraction close to 1/2 in the center of trap. For these 




FIG. 6: (color on-line) Top: The density for T/t = (left), 
T/t = 0.3 (middle), and T/t = 0.4 (right) for g/t = 1, u> = 
0.24, 6 P = 0, and 207 — 210 dipoles trapped. Bottom: x's are 
the checkerboard order parameter | (hi — hi+ ey ) | in the center 
of the trap as a function of T and o's are the checkerboard 
order parameter performed on the untrapped system at half- 
filling with the same parameters. 

parameters, there is a large region in the center of the 
trap with checkerboard density order for T = 0. With 
increasing temperature, the radius of the checkerboard 
phase in the center shrinks and it melts completely for 
T/t ~ 0.4. In Fig. [6] (bottom), we compare the cen- 
tral value of the density order parameter with that of an 
un-trapped system at half-filling performed on a 30 x 30 
lattice with the same interaction strength. We see that 
the critical temperature of the trapped system is close 
to that of an untrapped system. This shows that the 
system essentially behaves according to the local density 
approximat ion . 

In Fig. [7] (top) , we plot the density and the stripe or- 
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FIG. 7: (color on-line) Top: The density for T/t = (left), 
T/t = 0.1 (middle), and T/t = 0.3 (right) for g/t = 1, 
lu — 0.11, (9 P , 4> P ) = (tv/2, 0), and 179 - 190 dipoles trapped. 
Bottom: x's are the stripe order parameter \ {hi — fii+ ey }\ in 
the center of the trap as a function of T and o's are the stripe 
order parameter of an untrapped system at half filling with 
the same parameters. 



FIG. 8: (color on-line) The pairing order parameter 
(\{ci+e x &i)\ + |(c 4 _ ex c;}|)/2 in the centre for fi = —1.72 with 
6 P = tt/2, g/t = 0.85, and w = 0.11 as a function of T. There 
are 205 — 207 dipoles trapped. The inset shows a diagonal 
cross section of the pairing order parameter. Green o's are 
T/t = 0.05, red x's are T/t = 0.08, pink *'s are T/t = 0.09, 
and blue solid line is for T/t = 0.11. 



der parameter for the case where the dipoles are aligned 
along the x axis with (6p,<pp) — (tt/2,0). The coupling 
strength is g/t = 1, w = 0.11, and fi/t — —2 giving 
179— 190 dipoles trapped with the average filling / = 0.5 
in the center of trap. For this set of parameters, the cen- 
ter of the trap is in the striped phase for T = 0. The 
stripe order disappears with increasing temperature. In- 
terestingly, the melting is anisotropic in the sense that 
the stripe order disappears first in the y direction. The 
stripe order is completely gone for T/t ~ 0.3. Again, we 
see from Fig. [7] (bottom) that the density order in the 
center of the trap agrees well with that of an un-trapped 
system with the same parameters. 

Finally, we plot in Fig. [8] the pairing order parame- 
ter as a function of temperature for fi/t = —1.72 with 
(dp, (f> P ) = (tt/2, 0), g/t = 0.85, u = 0.11, and 205 - 207 
dipoles trapped. Since the coupling is weak, the system 
is superfluid for T = and there is no stripe order. As 
expected, the pairing decreases with increasing T and it 
disappears for T/t ~ 0.11. The critical temperature is 
calculated using mean-field theory. We expect correc- 
tions to mean-field theory to be small since the critical 
temperature is so small. The pairing increases slightly 
with increasing T at low temperature. This is because we 
for simplicity keep the chemical potential fixed in the nu- 
merical calculations leading to an increased density with 



increasing T. A number conserving calculation would 
yield a monotonically decreasing pairing with increasing 
T. 

These results illustrate that even in the presence of a 
trap, one can observe the superfluid and density ordered 
phases predicted for the infinite lattice systems, provided 
the system is large enough. In particular, the transition 
temperature is determined by the parameters in the cen- 
ter of the trap, and the results for a system with no trap 
can be used. 



CONCLUSION 

In conclusion, we examined the density ordered and su- 
perfluid phases of fermionic dipoles in a square 2D lattice. 
We determined the critical temperature of the density or- 
dered phases and demonstrated that it is proportional to 
the interaction strength for strong coupling. We calcu- 
lated the superfluid fraction and showed that the criti- 
cal temperature of the superfluid phase is proportional 
to the hopping matrix element for strong coupling. Fi- 
nally, we analyzed the effects of the harmonic trapping 
potential showing that for systems of a realistic size, the 
density ordered and superfluid phases exist with critical 
temperatures close to those obtained from a local density 
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approximation. 
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sions concerning the superfluid density and to S. Gam- 
melmark for Fig. [T] 
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